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ABSTRACT 

We suggest a physical mechanism whereby the acceleration time of cosmic rays 
by shock waves can be significantly reduced. This creates the possibility of particle 
acceleration beyond the knee energy at ~ \O i5 eV. The acceleration results from a 
nonlinear modification of the flow ahead of the shock supported by particles already 
accelerated to the knee momentum at p ~ p*. The particles gain energy by bouncing 
off converging magnetic irregularities frozen into the flow in the shock precursor and 
not so much by re-crossing the shock itself. The acceleration rate is thus determined 
by the gradient of the flow velocity and turns out to be formally independent of the 
particle mean free path (m.f.p.). The velocity gradient is, in turn, set by the knee- 
particles at p ~ p* as having the dominant contribution to the CR pressure. Since 
it is independent of the m.f.p., the acceleration rate of particles above the knee does 
not decrease with energy, unlike in the linear acceleration regime. The reason for the 
knee formation at p ~ is that particles with p > p* are effectively confined to the 
shock precursor only while they are within limited domains in the momentum space, 
while other particles fall into "loss-islands", similar to the "loss-cone" of magnetic 
traps. This structure of the momentum space is due to the character of the scattering 
magnetic irregularities. They are formed by a train of shock waves that naturally 
emerge from unstably growing and steepening magnetosonic waves or as a result of 
acoustic instability of the CR precursor. These losses steepen the spectrum above the 
knee, which also prevents the shock width from increasing with the maximum particle 
energy. 

Subject headings: acceleration of particles — cosmic rays — shock waves — supernova 
remnants — turbulence 



1. introduction 



The first order Fermi acceleration, also known as diffusive shock acceleration (DSA) is re- 
garded as the principal mechanism whereby galactic cosmic rays (CRs) are produced. The physics 
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of the energy gain is very simple: particles get kicked sequentially by scattering centers that are 
on the opposite side of the shock front and so approach each other. Its modern version has been 
developed by Krymsky (1977); Axford Leer and Skadron (1977); Bell (1978); Blandford and Os- 
triker (1978) and others. Soon after that, however, it was realized that the mechanism can only 
marginally (if at all, Lagage and Cesarsky 1983) account for the acceleration of the galactic cos- 
mic rays up to the knee energy at 10 15 eV (see also Jones et al. 1998; Kirk and Dendy 2001, for 
more recent good discussions). Up to the knee, the spectrum is nearly a perfect power-law so it 
is most probably produced by a single mechanism. At the same time, the low particle scattering 
rate makes the mechanism somewhat slow and the particle confinement to the shock insufficient to 
accomplish this task over the active life time of a typical galactic supernova remnant (SNR) shock, 
which are considered to be the most probable sites of CR acceleration. 

Not surprisingly, there have been many suggestions as to how to improve the performance of 
the DSA. In one way or another they targeted the most important and, at the same time, the most 
uncertain parameter that determines the acceleration time, the particle diffusivity K. The acceler- 
ation rate can be written roughly as %~} c = U^Jk{p), where U s h is the shock velocity. Therefore, 
the acceleration rate decreases with particle momentum p, since the particle diffusion coefficient 
K usually grows with it. In the Bohm limit, which is reasonably optimistic, one can represent K as 
K = cr g /3, where r g is the particle gyroradius (substituted in the last formula as a particle mean free 
path X) and c is the velocity of light. The condition X ~ r g is justified when the magnetic field per- 
turbations that scatter accelerated particles (and are driven by them at the same time) reach the level 
of the ambient magnetic field, 85 ~ Bq. Such a high fluctuation level has been long considered as 
a firm upper limit (McKenzie and Volk 1982) and there have been calculations indicating that the 
turbulence may saturate at a somewhat lower level due to nonlinear processes {e.g., Achterberg 
and Blandford 1986). 

It should be noted, however, that in cases of efficient particle acceleration, i.e., when a signif- 
icant part of the shock ram pressure pUf h is converted to the pressure of accelerated particles P c , 
there is indeed sufficient free energy of accelerated particles which can potentially be transformed 
to fluctuations with 85 ^> Bq. This would decrease X significantly below the particle gyroradius 
(calculated by the unperturbed field Bq), so the acceleration time would also decrease. From simple 
energetics principles, one can estimate the maximum fluctuation energy as (McKenzie and Volk 
1982) 



85 1 P c 

-2 ~M A -^ (1) 



where Ma = U s h/VA (Bq) is the shock Alfven Mach number. Bell and Lucek (2000) suggested 
that this upper bound may indeed be reached corroborating this idea by the numerical simulation 
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(Lucek and Bell 2000). The supporting simulation was rather limited in the system size and par- 
ticle energy range, which in turn, truncates the wave spectrum. Therefore, the saturation effects 
predicted earlier by McKenzie and Volk (1982); Achterberg and Blandford (1986); Blandford and 
Eichler (1987), have not been observed in the simulations by Lucek and Bell (2000). In any case, 
eq.(l) provides only an upper bound to the actual wave amplitude which may be much lower due 
to the saturation effects, not included in its derivation. 

A different approach to the same goal of increasing the turbulent component beyond the am- 
bient field level have been pursued in the recent papers by Ptuskin and Zirakashvili (2003, 2005). 
These authors studied the DSA in the presence of a Kolmogorov type MHD cascade initiated by 
the strong MHD instability in the CR precursor. Also recently, Bell (2004) took up from the Bell 
and Lucek (2000) and Lucek and Bell (2000) approaches and replaced the time consuming particle 
simulations by providing an estimate of the growth rate caused by the CR current at the shock. 
With the current kept fixed, the waves indeed grow to a high level sufficient to increase the max- 
imum particle momentum by a factor of 10 3 . It remains unclear, however, whether or not the 
CR trapping by the self-generated turbulence and other wave-particle interactions (e.g., Volk et 
al. 1984; Achterberg and Blandford 1986) saturates the instability at a lower level. Besides, the 
alleged growth of 85 beyond Bo clearly invalidates treatments of the linear instability which are 
based on the 85 <C Bq assumption. 

Diamond and Malkov (2004) suggested an alternative picture in which a strong 85 may be 
generated as a result of a type of inverse cascade in fc-space. Strong field perturbations driven 
resonantly by already accelerated particles are nonlinearly coupled to longer scales not only fa- 
cilitating the acceleration of higher energy particles but also providing an ambient field for them 
as the longest scale field perturbation. Such cascade can be driven by scattering of the Alfven 
waves on acoustic waves driven, in turn, by the Drury instability of the CR precursor. Another 
related route to long scale magnetic fluctuations is modulational instability of the Alfven waves 
themselves. The condensation of magnetic energy at long scales allows us to marginally preserve 
the weak turbulence 85 < Bq requirement, where the role of Bq is taken by the longest scale part of 
the magnetic energy spectrum (longer than the particle gyro-radius). This could possibly weaken 
the saturation factors setting in at 85 ~ Bo, which are noted above. 

The above-mentioned works, while aimed at the explanation of acceleration beyond the knee 
do not specifically address mechanisms that could be responsible for the formation of the knee 
itself. An exception is a paper by Drury et al. (2003), where the authors suggested that the knee 
can appear as a result of an abrupt slowing down of the acceleration at the sweep-up phase in 
combination with the Bell hypothesis about the generation of the strong magnetic field. Some 
further interesting ideas about the knee origin can be found in a recent paper by Sveshnikova 
(2003). 
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In this paper we propose a different scenario of faster than Bohm acceleration which is also 
intimately related to the knee phenomenon. Odd as it sounds, this mechanism does not require 
super-Bohm magnetic field fluctuations. The reason for that is that its rate does not depend on the 
particle scattering mean free path X. In the next section we describe the mechanism by comparing 
it with both the test particle and nonlinear DSA regimes under the Bohm diffusivity. We analyze 
the conditions under which the proposed acceleration regime is faster than the latter two. This 
sets the stage for Sec. 3 and 4 where a preliminary study of particle dynamics and transport in the 
CR precursor is presented. Section 5 deals with the estimates of the maximum energy that can be 
reached in excess of the knee energy. In Sec. 6 an acceleration model that allows us to calculate the 
slope of the spectrum between the knee and the maximum energy is presented. We conclude with 
a summary and brief discussion. 



2. Acceleration Mechanism: a primer 

Perhaps the easiest way to understand why and when the suggested version of the DSA be- 
comes faster than the standard one, is to consider why the latter is slow. For, we turn to the individ- 
ual particle treatment due to Bell (1978). Upon completing one acceleration cycle, i.e., crossing 
and re-crossing the discontinuity, a particle gains momentum 



p c 

where At/ is the relative velocity between the upstream and downstream scattering centers At/ = 
U\—U2~U\. Thus, over the acceleration time (when the momentum gain Ap ~ p), the number 
of cycles that the particle needs to make is N cyc i ~ c/U\ 3> 1. Apart from numerical factors and 
the differences between the upstream and downstream residence time contributions, the particle 
acceleration time can be estimated as 



x acc ~ ^ ~ Xc/Uf ~ x col c 2 /Uf (3) 

where X and x co \ are the particle mean free path (m.f.p.) and collision time, respectively. It is 
also useful to note that the acceleration time is of the order of time needed to the fluid element 
to cross the diffusion zone ahead of the shock, x acc ~ L^y (p) /U\, where L^/ = k(/>) j\J\ is the 
particle diffusion length. Using eqs.(2) and (3), one can write the following relation between the 
acceleration time, the time needed to complete one cycle and the collision time 

2 

C . C 

^acc '■ ^cycl '■ ^col ~ "j~fi * TT ' 
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Therefore, out of the c 2 /Iff wave-particle collisions, needed to gain a momentum Ap ~ p, only 
c/U\ are productive in terms of the energy gain. Most of the collisions are wasted. 

The situation changes fundamentally, when the number and energy of accelerated particles 
increase. First of all, the shock gets dressed with a cloud of accelerated particles (CRs) diffusing 
ahead of it, and the plasma flow becomes significantly modified by their backreaction. In addition 
to an abrupt velocity jump (which can be significantly reduced by this backreaction) the flow devel- 
ops an extended CR precursor (CRP, thereafter) of the length L p ~ k(/?*) /U\ = L^f (/>*), where 
p* is the particle momentum corresponding to the maximum contribution to the pressure of accel- 
erated particles. The flow in the CRP gradually slows down towards the main shock (subshock) 
and the spectrum is flatter than p~ 4 at high momenta, so that p* ~ p max and thus L p ~ L^f (/>*). 
Let us assume that a particle undergoes subsequent scattering by two scattering centers, approach- 
ing each other at a speed 8U <C c, Fig.l. Suppose that immediately after scattering off the right 
center, the particle has (in that center reference frame) the momentum p and the cosine of the pitch 
angle with respect to the plasma flow (moving to the left), fi < 0. After scattering off the left center, 
the particle momentum becomes p' and the cosine of the pitch angle // > 0, i.e., the particle moves 
back towards the first scattering center. Since scattering is elastic, from particle kinematics we 
obtain 



This result is obviously identical to the momentum gain in the conventional DSA theory (Bell 
1978). The left scattering center can be identified with one of the downstream, and the right 
one with one of the upstream, scattering centers. Considering /j and // as independent stochastic 
variables evenly distributed over the intervals (—1,0) and (0, 1), respectively, one can calculate an 
average momentum gain 

M. 2 % ) = ^ (5) 

p c 3 c 

where 

l r i 



Jo Jo 



Note, that the last result also holds in the case when there is an intervening scattering that does 
not reverse the particle direction. This can be seen from eq.(4) by noting that in such a case, both 
fi-s would be of the same sign and cancel upon averaging to the first order in 8U/c. Just as in the 
conventional (linear) acceleration theory a pair of scattering centers does not need (and is unlikely) 
to be the same at the next acceleration cycle. However, this is quite possible in the scattering 
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environment considered later in this paper. If so, then the momentum gain per cycle hp can be 
obtained from the conservation of the adiabatic invariant 



j> P\\dl\\ = const 



where the integral runs between the two scattering centers. If the interaction time of a particle with 
a scattering center is much shorter than the flight time between collisions, then 7|| ~ p\\l = const , 
where l(t) is the (decreasing) distance between the two centers. The formula (5) immediately 
follows from the last relation. 

The average time needed to complete the cycle is 

(&) = 2 /A\=4^ 

\Cfl/ c 

where we substituted X for an averaged distance between scattering centers. Now we can write 

dp_ (5p) _ 1 6U_ 

As long as X <C L p , one can approximate hU as 

dU 



8U ~ -X- 
oz 



so that the local acceleration rate is 



1 dU 

P^- 3 P^ (6) 

We observe that the particle m.f.p. dropped out of the acceleration rate in the smooth part of the 
shock structure, in contrast to the linear acceleration regime which occurrs at the shock disconti- 
nuity. The last formula is, of course, consistent with the standard diffusion convection equation 

d l + U d l^^ K ^l = X JJL p d l ( 7) 
dt dz dz dz 3 dz dp 

if one considers its characteristics ignoring diffusion effects. 1 Then, in the vicinity of the (abrupt) 
momentum cut-off p max , the main terms in eq.(7) are the first one on the left hand side and the 



'As it was shown by Malkov and Drury (2001), there is a reason to do so in certain situations. In particular under 
a strong nonlinear shock modification and Bohm diffusivity, (k « p) 7 the velocity profile is linear, i.e., dU /dz — const 
in a significant part of the CRP. 
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acceleration term on the right hand side. This clearly shows that a (front) solution propagates in 
momentum space at the speed given by eq.(6), basically independent of z and K. It is important 
to emphasize, however, that the independence of the acceleration rate of the particle diffusion 
coefficient (or m.f.p., as shown above) does not mean that the acceleration is faster than in the 
linear case where such dependence is the main factor making the acceleration relatively slow. 
The reason is very simple. The derivative dU /dz in eqs.(6) and (7) itself depends on K. Indeed 
estimating dU /dz as 

— — — — — (8) 

dz L p k(/?*) 

one sees that the acceleration time is x acc ~ k(/>*) ju\ i.e., particles with momenta p ~ p^ are 
accelerated approximately at the same rate as in the linear theory (apart from a numerical factor 
0(1), see Malkov and Drury 2001). Clearly, since the local (i.e., inside of the CRP) acceleration 
rate is independent of momentum, particles with p < p% accelerate more slowly than the linear 
theory acceleration rate, while particles with p > p* accelerate faster rate than in linear theory, 
and may in principle be accelerated faster than the Bohm rate. The overall acceleration process is 
illustrated in Fig. 2, depending on whether or not the pressure dominant momentum p* is stopped 
from growing after some t = t*. If it is, particle momentum grows exponentially, otherwise it 
continues to grow linearly in time. 

There are two problems with realization of this possibility. First, because in the nonlinear 
acceleration regime the spectrum at the upper cut-off is flatter than p~ 4 (it actually flattens to p~ 3 - 25 
at Pmax, Malkov 1997, in the limit of high Mach number and p max 3> mc), is clearly identifiable 
with Pmax- Therefore, unless the spectrum slope changes at p ~ p*, it is difficult to accelerate 
particles beyond p*. This is why the acceleration rates in both linear and nonlinear regimes are 
similarly slow, although for rather different reasons. In the linear regime, the acceleration rate 
slows down with time since the acceleration cycle duration, x cyc i grows with momentum (as gyro- 
period does), so the particle momentum grows only linearly with time. In the nonlinear regime, the 
width of the shock grows linearly with p*, so that the acceleration rate (8) decreases with p* (not 
with p\) and growth also linearly with time. Clearly, the first requirement would be to stop 
from growing at some point or at least to slow down its growth considerably. If = const, then 
particles with p > p* increase their momentum exponentially, eqs.(6,8). 

One simple reason for the saturation of (t) is a geometrical one. The thickness of the CR 
cloud ahead of the subshock, L p (p*) cannot exceed the accelerator size, i.e., some fraction of the 
shock radius R s , in a SNR, for example, (e.g., Drury 1983; Berezhko 1996). Thus we may identify 
p* with the maximum momentum achievable in SNR but limited by the geometrical constraints 
not by the lifetime of SNRs. What is also required here is that, along with or, independent of this 
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geometrical limit to p*, the character of particle confinement to the shock changes at when p* 
reaches some critical value. This should lead to a break in the spectrum at p ~ p* such that its 
slope at p > p* is steeper than p and flatter at p < p*. Then the main contribution to the particle 
pressure would come from particles around p*. 

For example, the spectral break may be initiated by a strong shock modification by accelerated 
particles (Malkov et al. 2002). Namely, Alfven wave compression in the CR precursor shifts 
the waves to the short wavelength end of the spectrum, leaving particles with p > />*without the 
resonant waves. In fact, some of those particles still remain resonant with the waves with k > 
r g l (P*) = ^* as l° n g as their pitch angle satisfy the condition p/j < m(i) c /k*, where p. is the cosine of 
the pitch angle, t0 c is the gyro-frequency. At the same time, the number of these resonant particles 
rapidly drops with momentum since their pitch angle distribution is limited by < mu> c /k*p, and 
they can hardly form a good pool for the acceleration. Rather, the momentum spectrum can be 
shown decay exponentially at p > p*. 

However, the latter example is based on a weak turbulence picture where a sharp resonance 
of particles with waves of random phases determines particle dynamics. As we demonstrate in 
the sequel, strong magnetic field perturbations result in a quite different character of particle con- 
finement. In the next section we consider a scattering environment which (though still simpli- 
fied) retains some important characteristics of a more realistic strongly unstable turbulent CRR 
As discussed earlier in this section, this turbulence should satisfy the following two conditions: 
i.) particles with momenta p > p* should not suffer catastrophic losses, ii.) their spectrum must 
nevertheless be steeper than p~ 4 . As we demonstrate in the next two sections, a gas of relatively 
weak shocks, that can emerge in a number of ways, provides the required scattering. 

3. Particle dynamics 

The conventional paradigm of particle transport in shock environment is the diffusive pitch 
angle scattering on self-generated, random phase Alfven waves. The resulting spatial transport 
is also diffusive with the coefficient that scales with the wave amplitude as K (p) bB k 2 , where 
the fluctuation wave number k must be set in the resonance relation with the particle momentum, 
kr g (p) ~ 1 . While this is the most consistent plasma physics approach that can be derived from 
the first principles using the quasi-linear theory provided that 85 <C Bq, there is no guarantee 
that it holds when the waves go nonlinear (e.g., Volk et al. 1984; Achterberg and Blandford 
1986). The latter is the rule rather than an exception in shock environments, which has been 
documented via observations, theory and numerical studies. So, the Earth's bow-shock as well 
as interplanetary and cometary shocks reveal a variety of coherent nonlinear magnetic structures. 
Usually, they originate upstream as a result of nonlinearly developed instability of the distribution 
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of ions reflected from the shock in the case of oblique propagation. If the shock is quasi-parallel 
the unstable ion distribution forms by thermal leakage from the downstream plasma. In cometary 
plasma, the pick-up ions drive the instability. Very often the magnetic structures are observed 
as an ensemble of discontinuities referred to as shock-trains or shocklets (Tsurutani et al 1986). 
There is no shortage of theoretical models describing these features, and we briefly return to them 
below. Standing somewhat outside of these models, but perhaps more specific to the nonlinear CR 
shocks is the acoustic instability driven by the gradient of the CR pressure in the CRR It was first 
studied by Drury (see e.g., Drury and Falle 1986; Zank et al. 1990; Kang, Jones and Ryu 1992, 
and references therein) and it also evolves into a gas of moderately strong shock waves propagating 
towards the main shock (sub-shock) in the CRR In what follows in this section, we will assume 
that the shock precursor is filled with the gas of such shocks separated by some distance L. 

Let us consider such a shock train propagating in the CRR Obviously the relative speed be- 
tween the individual shocks and the speed of the bulk flow is small compared to the particle speed, 
U (z) <C c. Therefore, we can consider the problem of particle scattering by scattering centers 
(shocks) and particle acceleration caused by the relative motion of the scattering centers, sepa- 
rately. Note that the acceleration part of the problem has been preliminarily considered in sec. 2 
and we will return to it in sec. 6. Turning here to the scattering part of the problem, we assume for 
simplicity that the shocks are one-dimensional and propagate at the same speed along the shock 
normal. Furthermore, we transform to their reference frame, neglecting compression by the main 
flow U (z) and their relative motion. Even under these simplifications, the shock train magnetic 
structures can be still rather complicated. For example, they have been extensively studied in the 
frame of the so-called derivative nonlinear Schrodinger equation and its modifications (Medvedev 
and Diamond 1995). Alfven - ion acoustic wave coupling results in wave steepening and gives 
rise to a shocktrain. A kinetic effect included in this model is particle trapping which was shown 
to be important for the formation of the shocktrain, which has been studied numerically (Hada 
et al. 1990; Medvedev et al. 1997). In another model, the shocktrain forms from a balance of a 
quasi-periodic driver and the nonlinearity of the magnetic field perturbations. The driver may orig- 
inate from one or a few unstable harmonics providing the smooth parts of the solution that steepen 
into shocks. These shocks also have a complicated spatial structure, typically characterized by a 
fast rotation of the magnetic field vector, caused by dispersion (see Malkov 1996, and references 
therein). However, since the gyroradius of high energy particles, even if reduced by possible mag- 
netic field amplification, is still much larger than both the width of the shock transition and the 
dispersive scale c/co p; , we replace each shock by a discontinuity of a coplanar magnetic field. A 
different possibility that results in a similar behaviour of the magnetic field occurs when the strong 
Drury instability discussed earlier, evolves into an ensemble of moderately strong shocks. 

Both cases are covered by the following representation of the magnetic field which corre- 
sponds to a periodic sequence of shocks, in which only one component of the field varies with the 
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coordinate in a co-moving reference frame, B = (0,B y , 1), where 

B y (z) =B y + B y sin [% (z - 1 /2)\ tanh | i cos [ji (z - 1 /2)] J (9) 

Here 5y is the constant component of the transverse to the shock normal component of the magnetic 
field and B y characterizes the strength of the shocks in the shocktrain (both normalized to the 
constant z-component, Bq), Fig. 3. We normalized the coordinate z to the distance between shocks, 
L. Of course, in reality different shocks in the shocktrain have not the same strength so that the 
coefficient B y should be replaced by a stochastic variable with a pdf, inferred from the shock 
dynamics. We return to this discussion in the next section, but in this simplified study we assume 
all the shocks in the shock train have the same strength. 

Let us consider particle dynamics in the magnetic field given by eq.(9). It is convenient to 
write the equations of motion using dimensionless variables in which time is normalized to L/c. 
Since we will be primarily concerned with the dynamics of particles having momenta p > p* 3> mc, 
we normalize their speed to the speed of light and set the absolute value of particle velocity to unity, 
V 1 . The remaining variable p is normalized naturally so as to the corresponding gyroradius 
r g (p) is measured in the units of L, that is the particle momentum p is scaled to eBoL/c. We also 
introduce the cosine of the pitch angle n = p z / p. The equations of motion read 



p = nxB (10) 
z = V (11) 

where n = p/p and p = const. 

Fig. 4 shows particle trajectories as a Poincare section, that represents orbit crossings of the 
shocks, i.e., the planes z = j (where j is an integer) on the gyrophase -pitch angle plane, (((),//). 
To make connection with the drift approximation we have transformed the particle momentum 
to the coordinate system rotated around the x axis by the angle f> = tan _1 fi y . In this coordinate 
system, particles spiral around the z axis with n = const if B y = 0. The gyroradius is taken to be 
equal to the distance between the shocks (p = 1, resonant case in terms of the linear theory) but 
the shock transitions are broad so that no distinct shocks with the steep magnetic field gradient 
are present and the particle dynamics remains largely regular. For smaller p the system becomes 
almost integrable since it can be reduced to a two dimensional one, because the particle magnetic 
moment is conserved with good accuracy. In the case shown in Fig. 4, however, some of the 
separatrices are clearly destroyed. Nevertheless, there is a clear separation of trapped (magnetic 
mirroring) and passing particles. Based on the conservation of the first adiabatic invariant / 
(l — fj 2 ) /B (magnetic moment), the trapped particles must be confined in momentum space to 
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Her — yj (B m ax ~~ Bmin) I B max ~ 0.27 which is close to what can be seen from Fig. 4. It is 
important to emphasize here that the trapping area does not shrink with p as was the case in the 
linear resonance situation, discussed at the end of the previous section. Clearly the spatial transport 
of these particles is non-diffusive. The trapped particles are convected with the shocktrain. The 
untrapped particles (|//| > /u cr ) propagate ballisiticaly and escape from the CR cloud. Thus, even a 
small amplitude (coherent) shock train provides a substantially different confinement of particles 
than the small amplitude random ensemble of Alfven waves does. It ensures perfect confinement 
of trapped particles and no confinement at all of untrapped particles. This is what is required for the 
acceleration mechanism outlined in Sec. 2. When this confinement regime takes over above some 
p = p*, since there are no resonant waves of the length larger than ^(p*), the energy spectrum 
must become steeper for p> p*. We will come back to this point in Sec. 6. 

The dynamics changes even more dramatically when the thickness of the shocks in the shock- 
train becomes small compared to the particle gyroradius. Taking the same parameters in eq.(9) 
as those used in Fig.4 except for v ~ 0, we show in Fig. 5 the Poincare section created by a sin- 
gle trajectory. It should be noted, however, that the actual amplitude of the shocks also increases 
somewhat in this case, even though the amplitude parameter B y is fixed. The phase space is now 
a "stochastic sea" with embedded "islands" of quasi-regular motion, quite typical for the deter- 
ministically chaotic systems (see e.g., Lichtenberg and Lieberman 1983). The dynamics inside 
the islands is weakly stochastic, so that the particle trajectories there are closer to the nearly in- 
tegrable case, shown in Fig.4, than in the rest of the phase plane. Therefore, within those islands 
where the averaged pitch angle Ji^O, particles propagate ballisticaly in the z direction, similar to 
the untrapped particles in the previous example. There are also islands where Jl « 0, and particle 
propagation is substantially suppressed there. The islands are stochastic attractors and the map- 
ping exhibits cycling around them apart from the quasi-periodic motion inside islands, not shown 
in Fig. 5. Note that only a single trajectory is shown. The first type of islands Qj ^ 0) are respon- 
sible for Levy flights (ballistic mode of particle propagation) whereas the second type of islands 
(Ji = 0) represent long rests or traps. The rest of the phase plane is covered by the region of a global 
stochasticity where particle propagation appear to resemble that around the islands but with shorter 
rests and Levy flights. Note that we use the term Levy flights for the above transport events even 
though the probability distribution function of the jump lengths might have a finite first moment 
(mean jump length). We investigate particle propagation in the z-direction in more detail in the 
next section. 
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4. Particle transport 

It is useful at this point to recall that in the standard diffusive shock acceleration the pitch- 
angle diffusion determines the spatial transport which is also diffusive. This picture results from 
the scattering of particles by Alfven waves of small amplitude and random phases. By contrast, 
in the preceding section we have considered a coherent magnetic structure of the finite amplitude 
shocks, as an alternative particle scattering agent at momenta higher than the pressure containing 
momentum scale p*. Despite the coherence, the particle dynamics, as we have seen, remains 
strongly chaotic and the onset of stochasticity depends on both the amplitude of the shocks and 
their thickness. We thus appeal here to the intrinsic stochasticity, caused by unpredictable particle 
motion in a field which can be perfectly regular. The fact that in a real situation the field is also 
random is often assumed to be of secondary importance. 

The resulting spatial transport is not a "classic" diffusion process. The connection between 
the particle dynamics, represented by the Poincare section shown in Fig. 5 and the particle spatial 
transport can be most easily understood by considering eq.(ll), z = i~i as a stochastic differen- 
tial equation, in which /j is a random process. Some general characteristics of this process can 
be inferred from the Poincare section shown in Fig. 5. As we mentioned above there are traps 
represented by islands which translate in the long flights or long rests depending on whether the 
averaged value of fu, i.e., Ji is nonzero or zero, respectively. Here the averaging is to be taken over 
an island attractor which is a layer of enhanced phase density around the island. 

The particle transport can be conveniently treated as a random walk on a lattice of shocks 
located at the integer values of z. However, in contrast to the classical random walk, this is a 
non-Markovian process since it is characterized by long rests (Zaslavsky 2002). Particles interact 
with the same shock repeatedly while they are trapped near an island with J2 = 0. Similarly, they 
perform a long jump in one direction and cross many shocks in a row while they are trapped in the 
phase space around an island with ju 7^ 0. 

The particle trajectory that corresponds to the Poincare section shown in Fig. 5, is demon- 
strated in Fig.6. The trajectory is indeed, non-diffusive. It consists of long 'Levy flights' that 
connect areas (clusters) of a rather slow particle propagation. Such clustered propagation regime is 
quite typical for nonlocal (fractional) diffusion models (see e.g., Metzler and Klafter 2000). The 
propagation within a typical cluster is magnified and shown in the inset of Fig.6. Certain similarity 
with the dynamics at large can be noted, which is also not unusual for the deterministically chaotic 
dynamical systems. However, the intra-cluster Levy flights are significantly shorter than the inter- 
cluster ones, the effect that will be shown to be much stronger for stronger shocks. Besides, long 
rests at different shocks can be seen inside the clusters. Note that it is this intra-cluster particle dy- 
namics that ensures their confinement to a localized region of the flow convected across the CRP 
as described in Sec. 2. As we argued there, every such passage through the CRP approximately 
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doubles the particle momentum. The Levy flights, on the other hand, can result in the escape of 
a particle from the shock precursor. We quantitatively describe and compare both possibilities in 
Sec.6.2. 

Let us come back to the overall dynamics. There are two important characteristics that de- 
scribe the random walk process. These are the waiting time and the jump length probability density 
functions (pdf). We have obtained these quantities from the trajectory shown in Fig. 6. The waiting 
time distribution at different shocks is shown in Fig. 7. Due to the gyro-motion of particles near 
the shocks, the pdf is not smooth but it has a smooth envelope which decays approximately as x~ 3 , 
where x is the waiting time. The distribution of the lengths of the Levy flights is presented in Fig. 8. 
Because of the asymmetry of the shock train (Fig. 3) there is a distinct, directional asymmetry of 
the jumps. Otherwise they exhibit an approximate power-law distribution up to the jumps hav- 
ing lengths limited by about 10. The pdf describes both the intra-cluster and inter-cluster particle 
jumps but only the jumps of the first type are suitable for the continuous acceleration. The longer 
jumps, as we estimate in sec. 5, generally result in the loss of particles. Another way of represent- 
ing the trajectory clustering is illustrated in Fig. 9, where the distribution of numbers of visits of 
different sites in the shocktrain is shown. 

The knowledge of the waiting time and jump distributions is necessary to derive the form 
of the adequate operators in the transport equation, (e.g., Metzler and Klafter 2000). Due to the 
non-Markovian character of transport these operators are integral (also called fractional differential 
operators, due to the memory in stochastic trajectories) rather than conventional diffusion opera- 
tors. In sec. 6 we shall take a somewhat simpler approach, which, however also requires the flight- 
rest time pdf to determine the transport coefficients. 

Up to now we have numerically examined the particle transport in shocktrains of a period 
longer or equal to the particles gyroradius. This corresponds to the resonant wave-particle inter- 
action in the small amplitude limit. Obviously even in the case p <C 1, there are resonant Fourier 
components in the shocktrain, since its spectrum decays as k~ 2 for k > 2n (recall that the shock- 
train period is set to unity). This is illustrated by Fig. 10, where particle trajectories with smaller 
momenta are shown. It should be also noted, however, that particles with smaller gyro-radii, such 
as P ^ P* can be efficiently confined due to the resonance with self-generated waves since, as we 
argued earlier, their population is significantly more abundant than those with p > p* because of 
the steeper spectrum in this momentum range. This result will be confirmed later. It is clear that, 
no matter how long the period of the shocktrain is, we need to estimate the confinement of particles 
with a gyroradius larger than that period, i.e., those with p > 1 . One can expect that their con- 
finement will be progressively deteriorated due to the fact that there will be only a high frequency 
force acting on these particles. To estimate the fraction of p 3> 1 particles that can be confined, we 
write eq.(10) in the form of the following system 
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p p\j\-[i L 

p = ^^Vl-^cos^) (13) 
p = const 

where we have assumed for simplicity that B y = 0. Because B y (z) has zero average, § and are 
particle gyro-phase and the pitch angle that are taken with respect to the shock normal. Assuming 
H 1 (which we verify below), from eqs. (11,13) we then obtain 

z-^B y (z)cos(<t>)=0. (14) 
Considering § as slowly varying compared to z and p, we can integrate the last equation once 

|z 2 +A (z) cos (4>) = const + <9 (p~ l ) (15) 

where we have introduced a periodic function A(z) according to B y (z) = —dA/dz. Since <f> is a 
slowly varying variable (it varies on a time scale p> 1, while z and vary on a time scale y/p), 
eq.(15) describes an oscillator lattice with a slowly varying potential as shown in Fig. 1 1 . The 
factor cos((|)) in eq.(15) slowly inverts the lattice potential. However, trapped particles do not 
detrap completely but rather get trapped in one of the nearby potential wells. The width of the 
trapping zone in z and thus in /j is given by \p\ < p s , where 




A m in) 1 



For By (z) given by eq.(9), we obtain: A min = —A max and A max = By/2%. Particles with > /li s are 
not confined to the shocktrain and propagate ballisticaly. 

To conclude this section, we emphasize that when the momentum of particles grows and 
their Larmor radius increases beyond the distance between the shocks, the fraction of particles 
that do not escape ballisticaly shrinks with momentum as p s <x \j ^fp, Fig. 12. This is slower than 
in the standard, weakly turbulent picture, where the critical value of p decays as I /p. On the 
other hand, the latter dependence is based on the linear cyclotron resonance without broadening, 
while the resonance broadening also improves the confinement (Achterberg 1981). In any case, 
one has to expect the spectral cut-off at p ~ 1 in units used in this section (Larmor radius is of 
the order of shock spacing). Therefore, the characteristic distance between shocks determines the 
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maximum momentum of accelerated particles. Of course, in a real situation there is a continuous 
strength/distance distribution of shocks which should produce a smoother spectrum decay at the 
cut-off momentum. Here we merely estimate the maximum distance which is equivalent to the 
minimum wave number in the standard acceleration picture. In the framework of the acceleration 
mechanism considered in this paper, the minimum wave number of the randomly phased Alfven 
waves should be of the order of r g 1 (p*), while particles with the momenta p* < p < 1 are confined 
via the mechanism of interaction with the shocktrain. In the next section we estimate the distance 
L between the shocks and thus p max since r g (p max ) — L, again, similarly to the standard estimate 
r g (pmax) — k min . 



In this section we estimate the maximum distance between the shocks in the shocktrain which, 
according to the preceding section, is directly related to the maximum momentum. In what follows 
we assume that the shocks are formed due to the development of an acoustic instability driven 
by the CR pressure gradient in the CRP. The latter has the scale height L p as being created by 
particles accelerated in a standard DSA manner and having maximum momentum p*. Specifically, 
L p ~ K(p*) jU\. As it was demonstrated by Drury and Falle (1986), the linear growth rate of this 
instability can be written as 



Here Pc and Pq z are the CR pressure and its gradient (which is antiparallel to the shock normal), 
respectively, is their adiabatic index, p is the plasma mass density and C s is the sound speed. 
For efficiently accelerating shocks, one can assume jc ~ 4/3. The first term in eq.(16) represents 
the wave damping caused by CR diffusion, calculated earlier by Ptuskin (1981). The second term 
is positive for the waves propagating along the pressure gradient while the oppositely propagat- 
ing waves are damped. The factor in the parentheses can reduce or even completely eliminate the 
instability, if dlnx/dp = — 1. However, there is no physical grounds, upon which this particular se- 
lection should be made (Drury and Falle 1986). It should be also noted here that particle diffusivity 
K refers here to particles with p < p* while the transport of higher energy particles is nondiffusive. 

Within the approximation leading to the growth rate given by eq.(16), there is no dependence 
on the wave number. A more thorough investigation can be found in Kang, Jones and Ryu (1992), 
which shows that the wave number dependence is indeed not strong. Under these circumstances 
it is reasonable to assume that the largest seed waves are the most important ones. Further we 
assume that the latter are related to the cyclotron instability of the Alfven and magnetosonic waves 



5. Estimate of the maximum momentum 




(16) 
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(which can provide the compressional seed component) that have excited by the accelerated parti- 
cles with momenta p < p.*, since particles of higher momenta have a steeper spectrum and do not 
contribute significantly to the growth rate of the cyclotron instability. We should focus primarily 
on the furthermost part of the precursor, since perturbations starting to grow there have the best 
chances to develop fully while convected with the flow towards the subshock. This part of the 
precursor is accessible only to the particles near the maximum momentum achievable within the 
standard acceleration scenario, i.e., p*. Hence, we assume that the typical wavelength of the Drury 
instability is Xo ~ r g (p*). 

While these perturbations grow to a nonlinear level, and propagate with the flow towards the 
subshock, they also steepen into the shocks. This behavior is seen in both two fluid simulations 
of Drury and Falle (1986) and, what is particularly relevant to our study, in the kinetic simulations 
of Kang, Jones and Ryu (1992). Due to the significant difficulties in numerical realization of the 
kinetic model, this study is limited in the maximum particle momentum, and thus in the length of 
the CRP, hp. Besides that, a monochromatic wave has been chosen as a seed for the instability. As a 
result, no significant interaction between the shocks has been observed. In a more realistic situation 
with a much longer precursor and shocks of varying amplitudes, strong interactions between shocks 
are to be expected (e.g., Gurbatov et al. 1991). In particular, stronger shocks overtake weaker ones 
and absorb them, so that the interaction has a character of inelastic collisions. Time asymptotically, 
such decaying shock turbulence is characterized by decreasing spatial density of the shocks. 

In the context of the present study, we must consider driven, rather than decaying, turbulence. 
A simple version of the driven Burgers model that reveals shocks merging and creation of new 
shocks has been considered by Malkov et al. (1995). In the case of the driven turbulence shocks 
not only merge but also new shocks are formed in-between. The new shocks merge with their 
neighbors and so on. In a steady state, some statistically stationary ensemble of shocks can be 
assumed with certain average distance L between the shocks and shock strength characterized by 
an average Mach number M. For a systematic study of particle propagation in such gas of shocks, 
the pdf 's of Mach numbers and distances between shocks are clearly required. As we mentioned 
above, we make here only rough estimates of these quantities. In any event, for the generation 
of statistically stationary shock ensemble, the instability should act faster than the convection of a 
fluid element across the precursor does, which requires JdL p /Ui 3> 1. The latter condition can be 
transformed to the following 

> 1 (17) 



pUiQ 

Note that the last condition does not contain the precursor scale L p and can be also represented as 
eM ^> 1, where £ = P c /piUf. With P c referring here to its subshock value, £ is the acceleration 
efficiency and M = U (z) /C s is an acoustic Mach number of the flow. Not surprisingly, the criterion 
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eM > 1 differs from the condition 85 2 / 5 o > 1 only by a factor of M/Ma = 1/a/P, as it follows 
from eq.(l). Indeed, the both instabilities share the same source of free energy, the CR pressure 
gradient. 

Now that we can assume that the Drury instability has sufficient time to fully develop to 
a nonlinear level we need to examine whether the nonlinear stage of the wave steepening lasts 
longer than the linear one. For the wave breaking time being shorter than the linear growth it is 
sufficient to fulfill Jo < Av/L, where Av is the characteristic magnitude of the velocity jump across 
the shocks and L is the characteristic distance between them. The last condition can be rewritten 
as 



±<^± (18) 
L p UieM 

Under these circumstances, in a steady state shock merging must be equilibrated with their growth 
due to the instability, i.e., y D ~ 8v/L, where 8v is the characteristic variation of the shock strength, 
i.e., the rms difference between the neighbor values of Av. The latter condition can be written as 



L Sv 1 

— ~ (19) 

L p UieM v 

Among the three parameters entering the last expression, dv/U\, £ and M, only the local Mach 
number M (z) needs special attention. Indeed, the far upstream M is considered to be prescribed 
and may be very large but because of the multiple shock formation in the CRP, strong heating 
should occur and drive this parameter down considerably. 

Let us estimate the local Mach number M. There are two major factors determining this 
parameter. After a fluid element passes through a shock in the shocktrain it gets heated but it 
cools adiabatically before it crosses the next shock. For this rough estimate we can assume that all 
shocks are of the same compression ratio and the overall flow pattern is L-periodic. For this order 
of magnitude type of estimate, we neglect the gradual variation of the flow at the scale L p > L. We 
thus assume that the flow speed in front of each shock is vi and it is equal to V2 = v\/r behind it, 
where r is the shock compression ratio. Similarly, the densities are related through p2 = rpi. We 
estimate the shock heating from the standard Rankine-Hugoniot relations. Let us denote the gas 
pressure in front of a shock in the shocktrain as P g \. Then, the gas pressure behind it is 

2 Y M 2 -y+1 
P * 2 ~ y+ 1 ^ 

The subsequent decompression phase decreases the pressure so that it becomes equal to P' gl = 
r~ y P g 2 in front of the next shock. Combining this with the above equation, we obtain 
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y-l+2/M 2 \ Y 2yM 2 -y+l 



Y+l J y+1 



Therefore the total pressure change after passing trough the shock and the smooth flow behind it 
is zero only for M = 1. However, the pressure does not really change for a broad range of Mach 
numbers around M = 1 (where the r.h.s. of the last expression is close to unity). This result means 
that M cannot be large, since it would lead to strong plasma heating driving its temperature to a 
state with M ~ 1 . However, in this situation the calculation of heating using the Rankine-Hugoniot 
relations is not correct, since shocks are likely to be sub-critical. Note that in our estimates, we 
also ignored the increase of magnetic energy behind the shocks. Since the above estimates show 
that shocks cannot be highly supercritical also, we conclude that a reasonable estimate for M to 
substitute in eq.(19) is a critical value of the Mach number i.e., M ~ M* which is of an order of a 
few (M* ~ 3, for example) depending on plasma (3 and f>„ g - the angle between the shock normal 
and the magnetic field (see Sagdeev 1966; Papadopoulos 1986, for a review of collisionless shock 
physics). 

The efficiency £ may then be calculated once the heating rate is known (see, e.g., Malkov and 
Drury 2001)). As for the parameter hv/U\, the first obvious constraint is 8v < Av which obviously 
verifies the consistency of eqs.(19) and (18). In fact, based on the studies of the driven Burgers 
turbulence (Gotoh 1998), we can assume that the shock strength variance 8v and its mean Av are 
related by 2 8v < Av. At the same time Av may be assumed to be not so much weaker than the 
subshock itself (Drury and Falle 1986; Kang, Jones and Ryu 1992). The subshock strength can 
be calculated from nonlinear acceleration theory (Malkov 1997; Blasi 2002; Blasi et al. 2005), 
along with the acceleration efficiency. Given the uncertainty of the heating rate one can estimate 
£~ 1/5 — 1/3, so that the estimate L/L p ~ 1 /5 does not seem to be totally unreasonable. 

As it follows from the previous section, the maximum particle momentum can be estimated 
from the relation r g (p max ) ~ L. Bearing in mind that L p ~ r g (p*)c/U\, the maximum momentum 
Pmax can be estimated as 



Since L/L p , as we argued, can be not very small and c/U\ is a large parameter, the suggested 
mechanism can produce a significant additional acceleration beyond the break momentum at p = 
p*. This stage of acceleration lasts for approximately x acc (p max ) ~ [ K (p*) /Uf] In (p max / P*)- 



Pmax 



c L 



UiL p 



2 This rough estimate should be taken with caution. Both quantities are really distributions, often with a power-law 
pdf, and may not even posses finite means. 
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6. Particle spectrum between the break and maximum momentum 

Up to now we have calculated the particle energy gain by their scattering on unspecified scat- 
tering centers carried by the converging flow in the CR precursor. Independently of that, we then 
considered the scattering of particles on an ensemble of shocks. The latter process did not lead 
to the energy gain, since no relative motion between the scatterers was assumed. Recall that we 
considered the simple shock ensemble as a magnetic structure traveling at a constant speed. Nev- 
ertheless, each shock in the shocktrain may posses its own flow structure with upstream and down- 
stream regions and thus can in principle accelerate particles via the standard DSA mechanism. 
However, as we pointed out earlier, these shocks do not have sufficient resonant turbulence up- 
stream and downstream with k < \/r g (p*). Therefore, there is no suffiecient coupling of particles 
with momenta p > p* to the converging flow around each such shock in order to gain momentum. 

In the situation we will consider further in this section, the momentum gain results from the 
combination of particle scattering off the scattering centers (shocks, localized magnetic structures) 
and gradual flow compression, leading to the convergence of the centers. In fact this is similar 
to what we have considered in Sec. 2 at an elementary level. An essential difference is that the 
scattering in momentum space is not homogeneous as the results of Sec. 3 suggest. Therefore, we 
relax the assumption, also made in Sec.2, that particle transport in momentum space is diffusion 
evenly covering isoenergetic surfaces. This affects the formula for the momentum gain. 



6.1. Particle momentum gain 

To calculate the rate at which particles gain momentum let us write down the equations of 
particle motion in the sub-shock reference frame and apply them in the area ahead of the subshock 
where the velocity of the flow is U(z,f) = U (z,t)e z (with z-axis pointing in the shock normal 
direction, e z ) 

? = -(v-U)xB (20) 
at c 

In this approximation, we have replaced the electric field by E = — c _1 U x B, assuming that due to 
the perfect conductivity the electric field vanishes in the local plasma frame. In contrast to Sec. 3, 
we do not specify the distance between the scattering centers here, so it is convenient to normalize 
the length to c/tO a -, time to CO" 1 , p to mc while v and U (z) to c. 

Introducing the cosine of the particle pitch angle to the shock normal, n, as well as the particle 
gyrophase (]), as in Sec. 4, it is useful to rewrite the equations of motion given by eq.(20) as 



-20- 



p = (be-*) (21) 

" = (22) 

* = — + -£^E=9t (23) 

z = ju ' (24) 

We have introduced the following complex variable instead of the transverse component of the 
magnetic field by b = (B x + iB y ) /B$. The above equations are written in the shock frame. Now 
we need to transform particle momentum p to the local plasma frame, i.e., the frame moving with 
respect to the shock with the velocity U (z, t) since it is this momentum the convection-diffusion 
equation [eq.(7)] refers to. This transformation can be written in the approximation U <C 1 as 
p' = (\—Up)p. Now we differentiate the both sides of the last formula with respect to time, 
taking into account that to this order of approximation dU /dt « zdU /3z. Making also use of 
eqs.(22) and (24), for the acceleration rate in the local plasma frame we obtain 

p' = ~» 2 -^p' (25) 

Recall that the precursor scale height L p that determines dU /dz should be calculated using the 
spectral break momentum at p*, i.e., L p = K(p*) jU\. We also note that the acceleration rate in 
eq.(25) does not contain the dynamical variable <|), and is actually almost independent of z. This is 
because U, to a good approximation, can be considered as linear function of z in a significant part of 
the shock precursor (in the case of the Bohm diffusion). As we mentioned earlier, due to the large 
gyroradius, we can neglect any short scale variation of U . This includes the structure of individual 
shocks in the shocktrain, so that only the gradual variation of U across the CRP is important. 
Also, we are really interested in averaged value of the acceleration rate p'/ p' . For example, if the 
particle transport in pitch angle is a small step diffusion, p can be regarded as random variable 
evenly distributed over the interval (—1,1). In this simple case, the average p? = j and we recover 
the standard acceleration rate given by eq.(6). In the case of structured phase space considered 
earlier in Sec. 3, one can assume that during most of an acceleration cycle, a particle is bound to 
the island having Jl = 0. If p ^ 0, the particle performs the Levy flights, which we discuss later. In 
the case of small islands (i.e., with the width <^ 1), we can replace p? by ju 2 (for Ji ^ 0) and by 
p 2 /3 (forju = 0). 

The situation becomes somewhat more complicated in the case of an oblique shock, where 
the phase space shown, for example in Fig. 3, implies that the cosine of the pitch angle p to be 
measured with respect to the average magnetic field and not with respect to the shock normal as in 



-21- 



eq.(25). The transformation to the angles ft and <j>, i.e., to the reference frame having z-axis aligned 
with the averaged magnetic field, is given by 



jj = jjcos-& — a/ 1 — ^sini^sin^j) 



(26) 



where f> = tan 1 B y (see Sec. 3). Therefore, the acceleration rate given by eq.(25) does include the 
gyrophase <j>. However, the dynamics on the plane in the case of trapping we are interested 
here and which is exemplified in Fig. 5, is mostly rotation around the origin, so that one can expect 
that Ji sin (j) ps 0. Thus the acceleration rate can be represented as 



P 
P 



-K 



„dU_ 

dz 



(27) 



where 



K 



^ocos 2 f3+^| 1 



1 



/jq sin f> 



(28) 



and /jq denotes the width of the island in which the particle orbit is trapped. We also omitted the 
primes in momentum p. 

Using this formula it is easy to calculate particle momentum gain assuming that the particle 
is trapped by the flow at the distance z from the subshock in the CRP where the local flow speed 
is U (z) and is convected with the flow to the subshock where the flow speed is Uq. Suppose that 
at the moment t and coordinate z the particle has the momentum pq. From eq.(27) we have for the 
particle momentum at z = where U = Uq 



p = po&xp 



K 



dt 



p exp 



if In 



U(z) 
Uq 



Po 



u_y 

Uq) 



(29) 



Note that in the conventional diffusive acceleration, where particles ergodically cover the entire 
isoenergetic surface, ie ji/o = 1, the acceleration constant K = ^ (independent of $) and the last re- 
sult signifies the effect of adiabatic compression of the CR gas by the converging flow. On the other 
hand, if the particle dynamics is such that /jq 0, $ ~ 7l/2, which corresponds to two-dimensional 
compression across magnetic field, one obtains K j. In strong CR modified shocks, the shock 
precompression R = U\/Uq scales as M 3 / 4 (Kazanas and Ellison 1986; Berezhko et al. 1996). As 

/ \4/3 

it was shown by Malkov (1997), the latter scaling is only valid for M < = ( V m jP*/p,„ 7 J 
and the shock pre-compression R saturates at R ~ Vi n jP*/Pinj when M > M* (see also Blasi et al. 
2005). Here the injection parameter Vi„j ~ (cpi„j/mUf) ncR/n\ and pt n j is the injection momen- 
tum while hcr and n\ are the number density of accelerated particles and that of the background 
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plasma far upstream, respectively. In both Mach number ranges the shock precompression can be 
quite significant unless the gas heating in the CRP is strong. At the same time, strong reduction of 
the shock precompression by plasma heating would diminish the heating itself by weakening the 
acoustic (Drury's) instability of the CRP. The situation may settle at some critical level where a 
moderate but still quite significant precompression is accompanied by the CRP heating caused by 
the acoustic instability (Malkov et al. 2000; Malkov and Drury 2001). Further interesting analysis 
of the nonlinear shock acceleration and its bifurcation has been published recently by Blasi et al. 
(2005). 



6.2. Acceleration model: details 



The picture that emerges from the above study of particle dynamics in the modified shock 
precursor can be described as follows. There are three groups of particles. One group is made up 
by particles that are locally trapped in the plasma flow by a structure in the shocktrain. They are 
convected with the flow, and either do not propagate at all or propagate very slowly with respect 
to the flow. They are clearly seen in e.g., Fig. 10 as the quasi-horizontal portions of the stochastic 
trajectory of a single particle. The remaining two groups of particles are particles performing Levy 
flights, i.e., those propagating ballisticaly in positive and negative direction at the speed U+ and 
U , respectively. These are characterized by steep portions of the trajectory in Fig. 10, with the 
positive and negative slopes, respectively. Let us denote the phase space density of the above three 
groups of particles by Fo(t,z,p) and F± (t,z,p). Furthermore, we introduce the probabilities of 
transition in the unit of time by denoting by a± the probability of trapping of particles performing 
Levy flight in positive and negative direction, respectively. We denote the rates of the reverse 
transitions by |3±. Now we can write the balance of particles in a phase space cell between zi, p\ 
and Z2 = z\ + dz, pi = pi+dp as follows 



^-F (t,z,p)dzdp = {F U\ Zi -FoU\ Z2 )dp+(pFo\ pi -pF \ p ^dz (30) 
- ($+ + $-)F dzdp + (a + F + + a-F_)dzdp (31) 

or after retaining only the linear terms in dz and dp and using eq.(27) we obtain the following 
equation for the phase space density of trapped particles Fq 

^ + luFo - K ^^P F o = ~ (P+ + P-) + a + F + + cc_F_ (32) 

Here the constant K is given by eq.(28). Proceeding in a similar manner we obtain two equations 
for F±, i.e., for the phase space density of the particles in the Levy flight state 
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dF ± d 



dU d 

K±^-^-pF ± = -a±F± + p±F 
az op 



(33) 



In contrast to eq.(32) the acceleration constants K± are defined by 



K± = 



\?± cos 2 ■& + - ( 1 - n\) sin 2 ■& 



(34) 



where p.± are the averaged values of the cosines of the particles pitch angles (jj± — fi) performing 
Levy flights in the positive and negative direction, respectively, while U± = c/j±cosf> 3> U are 
their speed components along the shock normal, cf. eq.(28). Note that we ignore direct transitions 
between the opposite Levy flights as comparatively rare events. This can be inferred from Fig. 10, 
for example. 

Equations (32) and (33) form a closed system that must be supplemented with the boundary 
conditions. Far upstream from the shock, at z = °°, it is natural to impose the following boundary 
conditions 



F (oo)=F_(oo) = (35) 

which simply mean that accelerated particles can only leave the system in that direction, but do not 
enter from it. Note, that F + (<*>) ^ and should be determined from eqs.(30) and (33). The second 
boundary condition, the one at the sub-shock location z = is somewhat more ambiguous. In the 
standard DSA scheme it is fixed by the flux of particles convected with the flow downstream. It 
can be expressed through the particle phase space density and the downstream flow speed since the 
particle distribution is isotropic. In the case of high momenta (p > p*) considered here, particles 
are not bound to the flow so strongly as to be able to make an isotropic distribution. Rather, accord- 
ing to eqs.(30) and (33), particles convected with the flow at the speed U (z) are being converted 
into particles propagating at high sped U± ~ c and can either leave the system upstream (unless 
they are re-trapped) or can penetrate deeply downstream and reach the contact discontinuity. The 
latter possibility for the high energy particles has been already discussed in the literature (Berezhko 
1996; Blondin and Ellison 2003). The fact that the contact discontinuity and the forward shock 
may indeed be closer to each other due to the backreaction of accelerated particles on the shock 
dynamics, as clearly follows from the nonlinear DSA, seems to find impressive observational con- 
firmation in (Warren et al. 2005). The strong magnetic field at the contact discontinuity may very 
well reflect particles with p > p* and they can return to the shock. In a simple form this require- 
ment can be obviously formulated mathematically setting the total particle flux through z = equal 
to zero: 
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F (0) U + F (0) U + F+ (0) U+ = (36) 

In other words, the negative part of the total particle flux, represented by the first two terms 
(Uo- < 0) is reflected by the strong magnetic turbulence downstream and ultimately reappears 
at the shock in the form of the untrapped particles propagating upstream at the speed U + > 0. The 
strong turbulence downstream may be formed by the magnetic perturbation convected and com- 
pressively amplified from the upstream media and by the Raleigh-Taylor instability of the contact 
discontinuity (June and Jones 1996). 

The way particle momentum enter eqs.(32) and (33) as well as the homogeneity of the bound- 
ary conditions given by eqs.(35) and (36), suggest the power-law solution: 

Fo,±<*p- q (37) 

Assuming for simplicity U± = const , (which is close to the truth, as e.g., Fig. 10 suggests) we can 
rewrite these equations as follows 



^-UFq = Yo^o + a+F + + a_F_ (38) 
dz 

dF+ 

U±—±- = Y±F ± + (3 ± F (39) 
dz 



Here we introduce the following notations 



Y± = -K ± (q-\) — -a± 

where 



P = P+ + P (40) 

From eqs.(38) and (36), we have 

Fo(z) = -^jfexp [u+F + (z') +a F (z')] dz' (41) 

To zeroth approximation in U/U± 1, from eq.(39) we have F + const and F const , while 
from the boundary condition given by eq.(35) and (36) we obtain 
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u + 

and F « 0. Upon substitution of F± into eq.(41) and taking it at z = we finally obtain the 
following equation for for the spectral index q: 

n^rw/o>)-<- 

In fact, the last formula allows a very simple physical interpretation, making it essentially equiv- 
alent to the well known general Fermi (1949) result. Indeed, Fermi showed that if a group of 
particles undergoes continuous acceleration at the rate x~} c and have certain probability to escape 
that can be characterized by the escape time x esc , then the energy (momentum distribution normal- 
ized above as Fq (p)) has the following spectral index 

q F = 1 + (43) 

tesc 

At least for some short period of time, by definition of these time scales, one can state that the 
momentum gain p(t)/p(0) = p/po = exp(?/x acc ). The probability to stay in the accelerator if 
escape is a Poisson process, t P C onf = exp(— t/x esc ). Then eq.(43) can be recast as (cf. Bell 1978; 
Blandford and Eichler 1987; Achterberg et al. 2001) 

In (l/'Pconf) 

In(PAPo) 

Note that t can be regarded here as the duration of an acceleration cycle and it does not enter the 
spectral index explicitly. It should be also clear that determination of Vconf and the momentum gain 
(or in other words, x esc and % acc ) requires in our model more information about particle transport, 
in particular, the particle trapping and untrapping rates a and p\ 

To demonstrate that the index q in our formula (42) must be interpreted in the same way as 
in the classical Fermi problem, let us analyze various factors in the integrand in eq.(42) separately. 
Assume that a particle enters an acceleration cycle at z = as a Levy flight in a positive direc- 
tion with a speed U+ against the plasma flow of speed U (z) in the shock precursor, located in 
a half-space z > 0. During its flight, the particle has a possibility to be trapped in the flow that 
is characterized by the rate a + . In other words, during time dt the probability to be trapped be- 
tween z and z + dz, equals oc + (t) dt = a + (z) dz/U + . We see this probability density in the integral 
given by eq.(42). Assume that the trapping event actually occurred at some z. Then the particle 
is convected (with the flow) back to the subshock at the speed U (z). During this convection it 
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has a certain probability to escape, which is characterized by the rate p\ This means that if the 
particle has a probability of being in the flow P(t), then the probability to remain there at t + dt 
is T (t + dt) = fP (t) (1 — $dt), or if the particle is certainly in the flow at t = (just trapped there, 
for example) then <£ (t) = exp (— Jq $dt) . Noting that dt = dz/U, we can write the probability of 
returning with the flow from z to z = as 



exp(-£$dz/U 



Combining with the above probability of being trapped between z and z + dz, we can express the 
probability of completing an acceleration cycle of the length between z and z + dz as 3 



dVcyd (z) = exp ^ jf jjdz^ oc+ (z) ^~ 



The total probability of completing one such cycle of any possible length < z < °° is 
We can now rewrite eq.(42) as follows 



Jo 



U(z) 
Uo 



cycl 



dz = 1 (45) 



dz 

Here the quantity G (z) = [U (z) /Uq] k is equal to a particle momentum gain p/fpi after it is con- 
verted with the flow from z to the origin, eq.(29). Using the mean value theorem, eq.(45) can be 
represented as 

if**- 1 

where (pf/pi) = G(z) is the mean momentum gain per cycle and P cyc i is, as stated above, the 
probability to complete this acceleration cycle. From the last formula we obviously have 

In ( 1 /Vcyci) 

<?=!+ : V 7 (46) 



3 Particle behavior during a Levy flight is non-Markovian, since it bears elements of deterministic motion (this is 
also true for the trapped particles). However, the joint probability of two such events is multiplicative since they evolve 
independent of each other. 



which is formally equivalent to the classical Fermi result. 



The quantities a + and B, which the index q depends upon, can be inferred from the particle 
stochastic trajectories shown, e.g., in Figs. 6 and 10. Indeed the probability for a particle to be in 
the trapped state can be written as % tr j (x tr + %l + x^ ) where the % tr and Xl are the average trapping 
and Levy flight times that can be calculated from the particle trajectory. Then, the particle trapping 
rate oc + can be written as 



(x ?r + x++x L )x+ 



similarly 



P± ~ 7 T =a 

[Xtr + H+ X L> % tr 

For the case of sufficiently strong shocks in the shock precursor, the trapping process is quite 
efficient, so we can assume that 



h B 

£</z«l. (47) 
o U 



Here the shock precursor length L p can be related to the trapping probability as 



L„~ / a + (z)zdz/ a + (z)dz. 
Jo Jo 

The requirement for the inequality in eq.(47) is then 



where the precursor crossing time is 



^ « 1 (48) 



[ l p dz 
Wo U 

An example of a particle trajectory with such long trapping times is shown in Fig. 13. Note that 
when a particle is bouncing between 2-3 neighboring shocks it can also be considered as "trapped" 
in the flow, as is the case for a particle stuck to a single shock. These kind of events are clearly 
seen in Fig. 13. 
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7. Summary and Conclusions 

The principal conclusion of this paper is that particle acceleration inside the cosmic ray shock 
precursor may very well be faster and proceed to higher energies during SNR shock evolution 
than the conventional DSA theory suggests. The requirement for such enhanced acceleration is 
some limitation on the momentum of particles that contribute the most to the CR pressure in order 
to prevent further inflation of the shock precursor. This can be achieved by spatial limitations 
of momentum growth, as discussed in detail by e.g., Berezhko (1996), by the wave compression 
and blue-shifting of the wave spectrum in the CRP (Malkov et al. 2002) and by the change of 
acceleration regime after the end of the free expansion phase (Drury et al. 2003). As a result of 
these limitations, a break on the particle spectrum is formed and maintained at p = p* beyond 
which particles do not contribute to the pressure significantly, but are still accelerated at the same 
rate as particles with momenta p < p*. 

The fundamental acceleration mechanism is essentially identical to the nonlinear DSA. The 

difference from the linear DSA is that frame switching which results in the energy gain occurs 
predominatingly between the scatterers convected with the gradually converging upstream flow 
in the precursor, and not between the upstream and downstream scattering centers, as usually 
assumed. The maximum momentum is estimated to be 



c L 

Pmax~ TT1-P* ( 49 ) 
U\ L p 

with the maximum value of L/L p < 1/5. Therefore, this mechanism provides up to ~ .2 (c/U\) en- 
hancement to the maximum momentum, based on the standard spatial constraints (e.g., Berezhko 
1996). The spectral index between p* and p max depends on both the ratio of the particle trapping 
to flight time % tr /%L m the upstream turbulent medium, and on the shock precompression R. This 
is in a deep contrast with the spectrum below which tends to have a form quasi-independent of 
R (e.g., Malkov and Drury 2001). The acceleration time is 



W- (Pmax) ~ ^nl (p* ) In (50) 

p* 

where x^l (p) — 4k (p) is the nonlinear acceleration time which only slightly differs from the 
upstream contribution to the standard linear acceleration time (Malkov and Drury 2001). 

The overall temporal development of the acceleration process can be described as follows. It 
starts as the standard DSA mechanism from some slightly suprathermal momentum pi n j and pro- 
ceeds at the rate x~ c \ (p max ) ~ Uf/K(p max ), where K(p) ~ cr g (p) and r g (p max ) k min ~ 1. The min- 
imum wave number k min decreases with growing p ma x(t), since the MHD waves confining particles 
with momenta p in j < p < p max to the shock are resonantly generated by the accelerated particles 
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themselves. After p max reaches a certain value (that can be calculated analytically, e.g., Malkov 
and Drury 2001) and is clearly seen in numerical time dependent solutions (e.g., Berezhko et 
al. 1996) a CRP forms and the acceleration then proceeds mostly in the CRP. The acceleration 
rate, however, behaves with time approximately as it does at the previous (linear) stage since 
dU jdz decreases with the growing width of the CRP, L p ~ K(p max ) /U\. Note that dU /dz sets 
the acceleration rate. The next acceleration boosting transition occurs when L p reaches one of its 
above mentioned "natural" limits, such as the accelerator size (a significant fraction of an SNR 
radius, for example, i.e., K.{pmax) /U\ ~ Rsnr)- Assume p max ~ p* at this point. Then particles 
with momenta p > p* cannot be confined to the shock diffusively, since their diffusion length 
Ldif(p) = K (p) /Ui > L p ~ k(/?*)/Ui. Fortunately, by this time instabilities of the CRP produce 
an ensemble of shocks or similar nonlinear magnetic structures. These ensembles of structures 
turn out to be capable of confining particles with p > p*, since the mean distance between the 
structures (which is the scale which replaces the wavelength of particle-confining turbulence in the 
quasi-linear picture) is much longer than the particle gyroradius, L^> r g (p* ) . Particle confinement 
is, however, selective in terms of particle location in phase space. This is similar to, but more 
general than, the particle confinement in magnetic traps where only particles having pitch angles 
satisfying \fi\ < ^AB/B are trapped while others leave the system freely. The price to be paid for 
this is a spectrum at p > p* which is even steeper than p~ A . This is in fact more a blessing than a 
curse since these particles do not contribute to the CR pressure significantly and L p does not grow 
maintaining thus the acceleration rate at the same level Uf/K.(p*) for all particles with momenta 

P*<P< p m ax- 

The total acceleration time up to the maximum momentum p max (given by eq.[50]) is thus 
only logarithmically larger than the acceleration time to reach the spectral break (knee) at p ~ p* 
and the maximum momentum itself is also pinned to the break momentum p* through eq.(49). 
The break momentum is simply the maximum momentum in the standard DSA scheme, whatever 
physical process stops it from growing further. The realization of the above scenario, which obvi- 
ously produces a significant acceleration enhancement in terms of both the maximum energy and 
acceleration time, depends crucially on the presence of the ensemble of shocks or similar scatter- 
ing magnetic structures in the CRP. Their existence in various (even CR unmodified) collisionless 
shocks is perhaps beyond any reasonable doubt, as we discussed in Sec. 3, but the real challenge is 
to calculate the statistical characteristics of their spacing L and amplitude parameter Av. A simple 
approach that we pursued in Sec. 5 leads to the conservative estimate of L/L p ~ 1/ 10, so that taking 
U\/c ~ 10 3 , the maximum momentum exceeds its standard value by a factor of 100, eq.(49). 

At this point it is worth while to recapitulate the major differences between our approach and 
other recent models discussed in the Introduction section which are also aimed at the enhancement 
of acceleration efficiency. They primarily seek to increase the turbulent magnetic energy by tapping 
the free energy of (already) accelerated particles. This would result in decreasing the particle 
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dwelling time (more frequent shock crossing, and thus faster acceleration) and, concomitantly, the 
particle diffusion length (better confinement, higher maximum energy). Our approach is based on 
the appreciation of 

(i) in a nonlinear regime, the acceleration time is reduced because of a steeper velocity profile 
in the smooth part of the shock transition (where acceleration mostly occurs). The steep gradient 
is maintained due to the formation of a break on the particle spectrum. Without the break, the 
gradient would flatten with growing particle energy slowing down the acceleration. Apart from 
this gradient, the acceleration rate does not depend on the m.f.p. explicitely, as long as the latter is 
smaller than the shock transition. 

(ii) particle confinement (maximum energy) is enhanced by increasing the wavelength (shock 
spacing in a shocktrain inside the CRP, for example) of turbulence which is in a distinctively strong 
regime. Departure from the weak turbulence is marked by modification of particle confinement in 
the CRP. This steepens the momentum spectrum beyond the spectral break and ensures its exis- 
tence, as required in (i). 

It remains unclear, whether the realization of the first (i.e., strong magnetic field amplification) 
hypothesis would preclude or otherwise strongly influence the acceleration scenario suggested in 
this paper. It would almost certainly not, if the magnetic energy increases predominantly at longest 
scales as suggested by Diamond and Malkov (2004). This should merely scale the estimates pro- 
vided in this paper by renormalising such basic quantities as L p and x NL , thus making acceleration 
improvement even stronger. One disadvantage of acceleration enhancement with significant field 
amplification, however, is that the increased magnetic field inevitably leads to increased radiative 
losses with a reduction of the particle maximum energy as an implication. This important accel- 
eration constraint has been studied in detail by Aharonian et al. (2002), including applications to 
many popular acceleration sites. Naturally, the most prominent impact is expected on the accelera- 
tion of the extremely high energy CRs (EHECR, > 10 20 eV). Detailed models of such acceleration 
have been proposed by Miniati (2001); Kang et al (2005) and others. Clearly, the acceleration im- 
provement without significant increase of the magnetic field potentially provides a significant edge 
on the field amplification scenario in the general context of the problem of EHECR acceleration. 

Support by NASA grant ATP03-0059-0034 and under U.S. DOE grant No. FG03-88ER53275 
is gratefully acknowledged 

REFERENCES 



Achterberg, A. 1981, A&A, 98, 161 



Achterberg, A., Blandford, R. D. 1986, MNRAS, 218, 551 

Achterberg, A., Gallant, Y.A., Kirk, J.G. and Guthmann, A.W. 2001, MNRAS, 328, 393 

Aharonian, FA., Belyanin, A. A., Derishev, E.V. Kocharovsky, V.V. Kocharovsky, VI. V, 2002, 
Phys. Rev., D66, 023005 

Axford, W. I., Leer, E. and Skadron, G. 1977 Proc. 15th Int. Cosmic Ray Conf. (Plovdiv) 11, 132 

Bell, A. R. 1978, MNRAS, 182, 147 

Bell, A.R. 2004, MNRAS, 353, 550 

Bell, A. R., and Lucek, S. G. 2000, Astroph. Space Set, 272, 255 
Berezhko, E.G., 1996, Astropart. Phys., 5, 367 

Berezhko, E. G., Yelshin, V, and Ksenofontov, L. 1996, Sov. Phys. JETP, 82, 1 

Blandford, R. D. and Ostriker, J. P., 1978, ApJ, 221, L29 

Blandford, R.D. & Eichler, D., 1987, Phys. Rep., 154, 1 

Blasi, P., 2002, Astropart. Phys., 16, 429 

Blasi, P., Gabici, S. and Vannoni, G. 2005, MNRAS, 361, 907 

Blondin, J.M., and Ellison, D.C. 2001, ApJ, 560, 244 

Chekhlov, A. and Yakhot, V, 1995, Phys. Rev. E, 51, R2739 

Diamond, PH., Malkov, M.A. 2004, J. Plasma Fusion Res. SERIES, 6, 28 

Drury, L. O'C, Aharonian, F. A., and Volk, H. J. 1994, A&A, 287, 959 

Drury, L. O'C., 1983, Rep. Progr. Phys., 46, 973 

Drury, L.O'C., van der Swaluw, E. and Carroll, O. 1983, astro-ph/0309820 

Drury, L. O'C. and Falle, S.A.E.G., 1986 MNRAS, 223, 353 

Fermi, E. 1949, Physical Review, 75, 1 169 

Gotoh, T. and Kraichnan, R.H. 1998, Phys. Fluids, 10, 2859 



Gurbatov, S.N., Malakhov, A.N., and Saichev, A.I. Nonlinear random waves and turbulence in 
nondispersive media: waves, rays, particles, Manchester University Press, Manchester 
(England) 

Hada, T., C. F. Kennel, Buti, B. and Mjolhus, E., J. 1990, Phys. Fluids, 11, 2581 

Jun, B.I., Jones, T.W. and Norman, M. L. 1996, ApJ, 468, L59 

Jones, T.W., Rudnick, L., Jun, B.I. et al. 1998, PASP, 1 10, 125 

Kang, H., Jones, T. W. and Ryu, D. 1992 ApJ, 385 193 

Kang, H., Jones, T. W. 2005, ApJ, 620, 44 

Kazanas, D., Ellison, D. C. 1986, ApJ, 304, 178 

Kirk, J. G., and Dendy, R. O. 2001, Journ. Phys. G., 27, 1589 

Krymsky, G. F. 1977, Dokl. Akad. Nauk SSSR, 234, 1306 (Engl. Transl. Sov. Phys.-Dokl. 23, 327) 
Lagage, P. O., and Cesarsky, C. J. 1983, A&A, 125, 249 

Lichtenberg, A.J. and Lieberman, M.A. 1983, Regular and Stochastic Motion, Springer (New 
York: Springer) 

Lucek, S. G. and Bell, A.R., 2001 MNRAS, 314, 65 

McKenzie, J. E. and Volk, H. J. 1982, A&A, 116, 191 

Metzler, R. and Klafter, J. 2000, Phys. Rep., 339, 1 

Malkov, M. A. 1997, ApJ, 485, 638 

Malkov, M. A., Diamond, P. H., and Volk, H. J. 2000 ApJ, 533, L171 
Malkov, M. A., and Drury L. O'C. 2001, Rep. Progr. Phys., 64, 429 
Malkov, M. A., Diamond, P. H., and Jones, T.W. 2002, ApJ, 571, 856 
Malkov, M.A., Physica D, 1996, 95, 62 

Malkov, M.A. , Kotelnikov, A.D. and Kennel, C.F. Physica D, 1995, 86, 480 
Medvedev, M. V., P. H. Diamond, 1995, Phys. Plasmas, , 3, 863 



-33- 



Medvedev, M. V. , Diamond, P. H., Shevchenko, V. I. and Galinsky, V. L. 1997, Phys. Rev. Lett., 
78, 4934 

Miniati, R, Ryu, D., Kang, H., and Jones, T. W. 2001, ApJ, 559, 59 

Papadopoulos, K. 1985 A Tutorial Review Geophys. Monogr. Ser. 34 ed. Stone, R.G. and Tsurutani, 
B.T. (Washington D C: AGU) p. 59 

Ptuskin, V.S. and Zirakashvili, V.N. 2003, A&A, 403, 1 

Ptuskin, V.S. and Zirakashvili, V.N. 2005 A&A, 429. 755 

Ptuskin, V.S. 1981 , Astroph. Space Set, , 76, 265 

Sveshnikova, I.G., 2003, A&A, 409, 799 

Sagdeev, R.Z., 1966 Cooperative phenomena and shock waves in collisionless plasmas, Rev. 
Plasma Phys. 4 (New York: Consult. Bur.) 

Tsurutani, B. T., R. M. Thorne, E. J. Smith, J. T. Gosling, and H. Matsumoto, J. of Geophysical 
Research, 92, 11074(1986). 

Volk H.J., Drury & McKenzie 1984, A&A, 130, 19 

Warren, J.S., Hughes, J. P., Badenes, C. et al. 2005, astro-ph/0507478, to appear in ApJ 
Zank, G. P., Axford, W. I. and McKenzie, J. F. 1990, A&A, 233, 275-284 
Zaslavsky, G.M. 2002, Phys. Rep., 371, 461 



This preprint was prepared with the AAS I^TgX macros v5.2. 



-34- 




Fig. 1. — Flow velocity profile in a CR dressed shock (solid curve). Hatched circles represent two 
scattering centers moving with the flow and approaching each other at a speed 8U. The filled circle 
depict a CR particle after the first and second collision (see text). 
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Fig. 2. — Schematic representation of acceleration process. Up to the time the acceleration 
proceeds at the standard rate, which is similar during both the linear and nonlinear phases (solid 
line). The particle pressure dominant momentum p* (t) grows as the maximum momentum of the 
entire spectrum. If for t > t* the pressure dominant momentum p* remains flat (solid line), particles 
with p > p* maintain the same acceleration rate that they had when their momentum was equal 
to /?* (dashed line). The dash-dotted line shows for comparison, how the maximum momentum 
Pmax = P* would continue to grow after t = U, had p* (t) not stopped growing at t = t*. 
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Fig. 4. — Particle orbits shown for v = 0.5; B y = 2; B y = 0.3 in eq.(9) and p = 1 on the (^/j plane. 
Each dot corresponds to the intersection of particle trajectories with one of the planes located at an 
integer z and coinciding with the shock locations (Poincare section). Several particle orbits form 
invariant curves (KAM-tori) as well as stochastic layers around some of the separatrices (see text). 
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Fig. 5. — The same as Fig. 4 except for the narrow shock transitions, v 0. In contrast to Fig.4, 
however, the entire phase portrait is formed by a single particle orbit, that does not visit the 'holes' 
in the phase space. 
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Fig. 6. — Particle trajectory represented as z (t) that corresponds to the phase space shown in Fig. 5. 
It starts close to the islands near /j = so that the orbit stays close to the origin (z = 0) for a long 
time. Such long rests will be repeated many times at different locations. Overall, particle transport 
is organized in clusters (five such clusters are shown) where the transport is strongly suppressed. 
The clusters are connected with long jumps where the transport is ballistic. 
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Fig. 7. — Waiting time distribution in a log-log format. Here n is the number of times the orbit has 
spent time % at any shock in the shocktrain. A somewhat discrete character of distribution is related 
to the gyromotion of particles near a shock so that when they interact with the shock repeatedly, 
the interaction time is commensurate with the gyro-period. 
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Fig. 8. — Distribution of Levy flights shown for forward and backward jumps. 




Fig. 9. — Distribution of number of visits of different shocks in the shocktrain. The magnified box 
shows the distribution inside of one of the clusters (see text). 
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Fig. 10. — The same as Fig. 6 but for smaller particle momenta, p = 0. 1 and p = 0.3. 
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Fig. 11. — Potential A (z) cos (<|>) in eq.(), corresponding to the magnetic field B y = —dA/dz, given 
by eq.(9) with v = and B y = 0. The function A (z) = (B y /n) [cos 71 (z - 1/2) - 1 jl\ . Solid line: 
cos ((])) = 1, dashed line: cos ((])) = — 1. 



-45- 




ln(p ) 



Fig. 12. — Phase space of particles in pitch- angle-momentum representation. The shaded area 
(Pinj < P < P*, \lA < 1) corresponds to the conventional particle confinement via randomly phased 
Alfven waves. For higher momenta, p > [if there were only weakly turbulent waves present 
with minimum wave number, k = k m i„ = \/r g (/?*)] this type of particle confinement would be 
limited to < /u c oc \/p, because of the resonance condition kr g (p)/j = 1. A shock-train with 
the spacing r g (p*) < L < r g (1) confines particles similar to the mirroring-type confinnement \n\ < 
He = a/AS/5 (hatched area), except for the phase space fragmentation (see Fig. 5 and text). Beyond 
p = 1 the particle confinement deteriorates to < /j s oc 1/ ^Jp~ (see text). 
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Fig. 13. — Stochastic particle trajectory represented in the same format as in Fig. 10 but for p = 0.3, 
By = 1 and By = 2. 



